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We numerically study the large deviation function of the total current, which is the 
sum of local currents over all bonds, for the symmetric and asymmetric simple exclu- 
sion processes with open boundary conditions. We estimate the generating function 
by calculating the largest eigenvalue of the modified transition matrix and by popu- 
lation Monte Carlo simulation. As a result, we find a number of interesting behaviors 
not observed in the exactly solvable cases studied previously as follows. The even and 
odd parts of the generating function show different system-size dependences. Different 
definitions of the current lead to the same generating function in small systems. The 
use of the total current is important in the Monte Carlo estimation. Moreover, a cusp 
appears in the large deviation function for the asymmetric simple exclusion process. 
We also discuss the convergence property of the population Monte Carlo simulation 
and find that in a certain parameter region, the convergence is very slow and the gap 
between the largest and second largest eigenvalues of the modified transition matrix 
rapidly tends to vanish with the system size. 
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1. Introduction 

Current, such as electric current, heat flow, or current of matter, plays important 
roles in nonequilibrium systems. The presence of macroscopic current is a signature of 
nonequilibrium states. Near equilibrium, the current is proportional to the conjugate 
external field and the linear coefficient is given by the autocorrelation function of the 
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current in equilibrium states. Beyond the linear regime, although we do not know 
much about general properties of the current in strongly interacting systems, current 
fluctuations are considered as key quantities to understand the dynamical behavior of 
nonequilibrium systems. In this respect, the large deviation property of the current 
is being studied intensively^"^-* in various systems. The large deviation function char- 
acterizes the probability of the time-integrated current in the long-time limit and its 
Legendre transform gives the cumulant generating function. Thus, it contains informa- 
tion about not only the mean and the variance but also higher-order cumulants, which 
are related to the nonhnear transport coefficients in the full counting statistics.^) In the 
context of statistical physics, the time-integrated current called the Helfand moment 
is related to the Green-Kubo formula.^) The shear viscosity and the thermal conduc- 
tivity calculated with the Helfand moment in the recent simulations are in excellent 
agreement with those calculated with the Green-Kubo formula. ^'^^ It is also suggested 
that the current large deviation is relevant to understanding the dynamical behavior of 
glasses. 

The simple exclusion process (SEP) is probably the simplest and best studied model 
of nonequihbrium systems. ^°'^^) The SEP is defined on the one-dimensional lattice, 
where particles are allowed to move to one of the nearest neighbor sites when the site 
to move to is empty. If the leftward and rightward hopping rates are the same, the 
model is called the symmetric simple exclusion process (SSEP). If the hopping rates 
are different, the model is called the asymmetric simple exclusion process (ASEP). Two 
extreme cases of the ASEP have special names. When the difference is as small as the 
inverse of the system size, the model is referred to as the weakly ASEP (WASEP). The 
totally ASEP (TASEP) means that one of the hopping rates becomes zero. The ASEP is 
applied to various phenomena such as traffic ffow,^^-* molecular motor transportation,^^) 
the process of copying messenger RNA,^^) and sequence alignment. The ASEP in 
the infinite system with a step initial condition can be mapped to a model of surface 
growth, and recent experiments show that the height fiuctuations of the growing surface 
in electroconvection, which correspond to the current fiuctuations in the ASEP, obey 
the Tracy- Widom distribution function that appear in the exact solution of the one- 
dimensional Karar-Parisi- Zhang equation. 

For the ASEP with open boundary conditions, the stationary state is exactly cal- 
culated with the matrix product method^^"^^-' or with the recursion relation. In par- 
ticular, it is clarified that the system has three phases, the low-density phase, the 
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high-density phase, and the maximum-current phase, depending on the particle's enter- 
ing and exiting rates at the boundaries. The transition between the maximum-current 
phase and the other phases is second order and that between the low-density and high- 
density phases is first order. Thus, in the latter case, the two phases can coexist on 
the phase transition line, where a kink appears between the low-density region and the 
high-density region. The random walk motion of the kink contributes to the power-law 
behavior of the power spectrum. ^^'^^^ 

The current large deviation functions for the SSEP with the open boundary condi- 
tion were studied by Derrida et al.,^'^^^ who obtained the scahng form for the generating 
function. The TASEP^^'Se) ^nd the WASEP^^) on a ring, and the SSEP in the infinite 
system with the step initial condition^^^ were studied by applying the Bethe ansatz. 
The direct calculation of the probability of the number of particles passing through 
a certain site was performed on the TASEP and the WASEP in the infinite system 
with some special initial conditions, and it turned out that different initial con- 
ditions lead to different types of scaling behavior. The current moments of the ASEP 
in the infinitely large system have recently been studied by Imamura and Sasamoto,^^^ 
who showed that the calculation of the nth moment is reduced to the problem of the 
ASEP with a system with particles less than or equal to n. The macroscopic fluctua- 
tion theory has been developed as a theory for the large deviation function. ^^"^^^ It is 
based on the postulate that the large deviation function is characterized only by the 
mean and variance of the current. The macroscopic fluctuation theory is applied to 
the SSEP^®) and the WASEP. In the ASEP, the joint probability of the current and 
the density is calculated for some restricted boundary conditions by using the matrix 
product method. For the current through the boundary, the generating function in 
the ASEP with open boundary conditions where both input and output are allowed at 
each boundary is calculated by using the Bethe ansatz^^^ for the low- and high- density 
phases, and the symmetry relation for the current is also discussed. 

In addition to the exact calculations, many numerical studies have been carried 
out. A population Monte Carlo method^°) was devised to calculate the large deviation 
functions flrst for discrete time dynamics,^^^ and it was extended to continuous time.^^^ 
By using these methods, the large deviation functions in the zero range process^^^ 
and in the Kipnis-Marchioro-Prcssuti (KMP) model, '^"^^ which is a simple model of 
heat conduction, were computed. The density matrix renormalization group (DMRG) 
method was also applied to obtain the current large deviation for the TASEP with open 
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boundary conditions.^^^ 

Despite these extensive studies and the importance of the apphcations, the current 
large deviation function in the ASEP with open boundary conditions has not yet been 
fully understood. Furthermore, the definition of the current varies in the studies. Some 
studies consider the current through a bond and other studies consider the total current, 
which is the sum of the former over all bonds. Although the averages of these two 
quantities are trivially equivalent, their fluctuations may behave differently. 

In the present study, we numerically calculate the current large deviation function 
in the SEP with open boundary conditions and find several properties that have not 
been reported in the exactly solvable cases. The large deviation function shows a cusp 
near zero current when the asymmetry is large. In both the SSEP and the ASEP, the 
even part of the generating function depends linearly on the system size, while the odd 
part does not. Our finding indicating that the even-order cumulants are proportional to 
system size L differs from that obtained in the preceding study, where the nth-order 
cumulant is proportional to L"^^. We compare the generating functions for the two 
definitions of the current and find that they coincide with each other in the direct eval- 
uation using the largest eigenvalue method, while the population Monte Carlo method 
fails when the current through a boundary is employed. Thus, the use of the total cur- 
rent has an advantage in the Monte Carlo simulation. We estimate numerical errors in 
the Monte Carlo simulation using the symmetry relation and find that the deviation is 
usually small. However, the Monte Carlo simulation sometimes shows an unexpectedly 
slow convergence and produces no reliable results. In those cases, the difference between 
the largest and second largest eigenvalues of the modified Master equation decays faster 
than the exponential of the system size. 

This paper is organized as follows. In §2, we give a brief review for the SEP and the 
current large deviation function, and explain two numerical methods that we use in the 
present simulations involving the use of the largest eigenvalue of the modified Master 
equation and the population Monte Carlo method. We also give a short summary of the 
studies of the current large deviation for the SSEP and the ASEP under open boundary 
conditions. In the next section, we show our numerical results. In particular, we focus 
on the difference induced by different definitions of current. The convergence problem 
in the population Monte Carlo method is also discussed. The last section is devoted to 
discussion and conclusions. 
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2. Model and Algorithm 

2.1 ASEP and Master equation 

The SEP is a continuous-time Markov process defined on a one-dimensional chain. 
We denote Tj = 1 if site j is occupied by a particle and Tj = if the site j is empty. 
The configuration of the system of size L is described by the set C = {ti,T2, . . . ,tl}. 
A particle can hop to the left or right nearest-neighbor site if the destination site is 
empty. Namely, the particles are subject to the exclusion interaction. The hopping rate 
to the right is denoted by Pr and that to the left by pi. If site 1 is empty, a particle 
enters with rate a, and if the site is occupied, the particle exits with rate 7. Similarly, 
if site L is empty, a particle enters with rate 5, and if the site is occupied, the particle 
exits with rate f3. 

Let us denote the probability for the system to be in configuration C at time t by 
P{C,t). By using the transition rate from the configuration C = {r[, . . . ,r|} to C = 
{ti, . . . , tl}, Wcc'-i the time evolution of P(C, t) is described by the Master equation 

^P(C,i)= [Wcc'P{C\t)-Wc'cP{C,t)]. (1) 
In the ASEP, the transition rate is written as 

WcC = WiSr^r^ ■ ■ ■ 6r^r'^^ 

L-l 

j = l 

+<^rir{ •••(^ri_iri_,'ti'L, (2) 

where 5^^^/ = (1 - Ti){l - r/) + r^r/, 

wi = a(l - t[)ti + ^t[(1 - Tl), (3) 

WL = 5{l-r'L)TL + PTi{l-TL), (4) 

and 

Wjj+i = PrTj{l - Tj+i)(l - Tj)Tj+i 

+Pl{l-T'X+^T,{l-T,+^). (5) 

If we denote r(C) = ^C'C, eq. (1) is rewritten as 

^P(C,i)= J2 Wcc'P{C',t)-r{C)P{C,t). (6) 

Thus, r{C) means the rate of transition from C to any other configurations. 
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The time evolution of the system is decomposed into the set of one-particle move- 
ments that we call a path. Consider that the system starts with initial configuration 
Co at time to. The configuration changes from Cj_i to Cj when a particle changes its 
position or a particle enters or exits from the system at time tj {i = 1, . . . , K), reaches 
Ck at tK, and remains in the same state until time t. Thus, the state at time t is 
C = Ck- This path is represented as (C, t; C^, ti^; . . . ; Ci, ti|Co, to)- Now, we denote 
the probability that the transition from Cj_i to Cj occurs between time tj and tj -|- dtj 
for i = 1, 2, . . . , X by 

P(C, t; Ck, tK', ■ ■ ■', Ci, ti|Co, to)dtA- . . . dt2dti. (7) 

Then the path probability density P{C, t; Ck, tK', ■ ■ ', Ci, ti|Co, to) ior K > 1 is written 
as 

P(C, t; Ck, tK', ■ ■ ■ ', Ci, ti|Co, to) 

. . .g-r-(Ci)(t2-ti)|y^^^^g-r(Co)(ti-to) 



Scck exp 



K 



K-l 



U^c,,.o„ (8) 

j=0 



-J2r{Ci){ti+i-ti) 

_ i=0 

where tK+i should be understood as t. Because K represents the number of one-particle 
movements, if K — 0, P{C, t; Ck, tK', ■ ■ ', Ci, ti|Co, to) should be interpreted as repre- 
senting the probability that no transitions occur at the time interval [to, t], which equals 
ScCo^~^^'^°^^^~^°^ ■ The relationship between the conditional probability P(C, t|Co, to) and 
the path probability density is given by 

P(C,t|Co,to) 

^ J2 Yl / dtA / dti^_i • • • / dti 

K=0Ci,C2,...,Ck^° *° *° 

P(C, t; Ck, tK', ■ ■ ', Ci, ti|Co, to). (9) 
If we differentiate the above equation with respect to t, we reobtain the Master equation. 



2.2 Large deviation function and generating function 

Let us denote the probabihty that the mean current takes the value q as 



K-l 

q^-J2QiCi+i,C,) (10) 



1=0 
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by -P(g). The microscopic current Q{Ci+i,Ci) takes the value 1(— 1) when the particle 
moves forward (backward) during the configuration change from Cj to Cj+i. A precise 
definition of the microscopic current is given in the next subsection. When t is large, 
this probability asymptotically behaves as 

P(g)~e-*^(«\ (11) 

where /(g) is called the large deviation function. As a result of the law of large number, 

f{q) must have a minimum at a certain value q* and satisfy f{q*) = 0. The Legendre 
transform of f{q) is called the generating function, which we write 

/x(A)=max[Ag-/(g)]. (12) 

The generating function satisfies fi{0) = that corresponds to f{q*) = 0. Similarly, the 
large deviation function is given by the Legendre transform of the generating function 
as 

/(g)=max[gA-/.(A)]. (13) 

A 

The generating function is also rewritten as 

for large t. This means that the generating function is expanded as 

10. lb. 

n=l n=l 

with respect to the cumulants {q'^)c, or the cumulants are generated as 

1 9"/x(A) 



(16) 



A=0 



which is the origin of the name of the function. 

The expectation in the rhs of eq. (14) is taken with respect to the path probability. 
Thus, we obtain 



= E E / dtic / dtK-i--- / dti 

§^,^^^^T.f=o^QiCi+uCi)^-riCK)it-tK) 
K-l 

n W^c.+iC,e-^(^*)(**+^-*^). (17) 



i=0 

If we introduce W^(j, = Wcc^-'^^''^''^"' and 



p\c,t\CoM 
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K=nnr^ ... r7i •^*o J to J to 



K=OCk,-,Ci 

5c.ce^^-^(^')(*--*') n (18) 

i=o 

the generating function is obtained as 



//(A) = lim - log 

t—^oo t 



J2P\C,t\Co,to) 



(19) 



c 

or we may assume some initial distribution Po(C'o) and substitute 

P\C, t) = J2 P\C, t|Co, to)Po(Co) (20) 

Co 

in place of P'^(C, i|Co, to) in the rhs of eq. (19). The time evolution of P^{C,t) obeys 
the modified master equation 

^P\C,t) = W^c'P\C',t)-r{C)P\C,t). (21) 

We note here that because '^c'i^c) ^cc "''(C*); this equation cannot be interpreted 
as an evolution of probability. Namely, X^c" -P^(C', t) varies in time. In the vector form, 
eq. (21) is written as 

ApA(i)^W^P(t), (22) 

where P^(t) is the vector whose Cth element is P^{C,t), and the A-modified transition 
matrix is defined as 

W = < ■ 23 

^ '^^ I -r(C) for C = C" 

Because eq. (21) is linear, the solution is spectrally decomposed into 

P\C^) ^J2^i-'MC)J2'^'n{C')Po{C'), (24) 

n C 

where (n is the eigenvalue and ip'niC) and ipn{C) are the corresponding left and right 
eigenvectors of W'^. In the long time limit, only the term with the largest eigenvalue Cmajx 
is dominant, and accordingly the generating function becomes equal to Cmax- Thus, we 
can estimate the generating function by numerically calculating the largest eigenvalue 
of the A-modified transition matrix if exact calculations are not possible. This method 
is powerful but restricted to small systems because of memory limitations. 
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2.3 Definition of the microscopic current 

There are two different definitions of current. One is the current through the bound- 
ary of the system. From the equation of continuity, this current determines changes of 
the amount of the conserved quantity in the system. The other is the total current given 
as the spatial integration of the current density. It is natural to use the total current 
in the linear response theory when we consider the response to a static and uniform 
external field or thermal force. 

Both definitions appear in the hterature of SEP. In the studies of current large 
deviations of the SEP with periodic boundaries, the total current is employed by Derrida 
and Lebowitz,^^) Derrida and Appert,^^^ and Prolhac and MaUick.^^^ On the contrary, 
the current through the left boundary is used by Derrida et al.^'^^^ or de Gier and 
Essler,^^'^^^ who study the systems with open boundaries. 

We let Qa denote the total current defined as 



L-l 



- a - rr)r^' + rlil - r^^') 
+(i-ri)rt'-rl(l-Tt'), (25) 
and Qb denote the current through the left boundary as 

Qb{Q+,,Q) = -(i-rl)rl+' + Tl{l-Tl+'). (26) 
For each current, we define the mean current qx {X — A or B) as 

K-l 

dx^-Y^Q^i^i+^^^i)^ (27) 

and the generating function //^(A) as 

gMx(A) ^ {e,MxtY (28) 

We note that the positive direction of the currents (25) and (26) is defined as leftward 
in this paper. We use the character q for a current if we do not need to specify qA or g^, 
or for a current in general. Since the system goes to a stationary state in the long-time 
limit, we obtain 

{qA) = {L + l){qB). (29) 
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Thus, the two definitions are equivalent with respect to expectation values except the 
trivial factor. 

Our question here is whether the large deviation properties of the two currents 
appear the same or not. Theoretically, they arc considered equivalent as follows. '^^^ If a 
path is fixed, most particles that enter the system through a boundary exit from either 
boundary and the contribution from such particles to is exactly L + 1 times the 
contribution to Qb- The number of such particles is proportional to t. On the contrary, 
the particles initially present in the system or those remaining in the system at time t 
contribute to and qB differently. However, the number of those particles should be 
^-independent and the contribution from those particles is negligible for a large time. 
Thus, qA and qB are supposed to be essentially the same. Despite such consideration, 
how they appear in the numerical experiment is a different problem. Actually, as will 
be clarified later, the use of the total current has an advantage in the Monte Carlo 
simulations. 

2.4 Monte Carlo simulations 

If the system size is so large that the direct evaluation of the largest eigenvalue 
is not possible, we can employ a Monte Carlo method. However, because Px{C,t) is 
not a probability distribution function, conventional Monte Carlo methods cannot be 
used. In fact, an algorithm based on the diffusion Monte Carlo method is devised by 
Giardina et al.^^^ and extended to a continuous-time case by Lecomte and Tailleur.^^^ 
In those algorithms, besides the conventional Monte Carlo dynamics, we introduce 
clones of configuration C. The clones are duplicated or pruned with the rate y{C) — 
Q{rx{C)-r{C))At ^ whcrc At is the Poisson-distributed waiting time and its distribution 
function is given as 

rA(C)e-'-^(^)^*, (30) 

where 

rx{C) = J2 ^cc- (31) 

Thus, the mean waiting time is l/rx{C). Since the number of clones must be limited in 
the computer, we actually need an additional algorithm to keep the number of clones 
constant. If the clone is to be pruned, we choose a random clone and copy over the 
pruned one, and if the clone is to be duplicated, we copy the clone on the randomly 
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chosen clone. 

The precise algorithm is as follows. 

(1) Set initial conditions for N clones. 

(2) Choose which clone to evolve. Here, we name it c". First, the clones are chosen 
orderly. After all the clones evolve in the initial step, the clone with the earliest 
time is chosen. 

(3) Calculate the transition probability W^(.,/rx{C), and let a transition take place. 
The time interval At is chosen from the Poisson law. 

(4) Calculate y{C) and set y = [y{C) + ^], where ^ is a uniform random number in 
< ^ < 1 and [x] is the integer part of x. If y = 1, the clone is preserved; if 
y — 0, the clone is erased and overwritten by another randomly chosen clone; 
and if y > 1, y — 1 clones are chosen randomly and overwritten with the clone c". 

Thus, the total number of clones is kept constant. We denote y at the ith time step 
by yi and define Xi = {N + yi — l)/N. Thus, the generating function is given by the 
following formula in the long-time limit i^^-* 

li{X)^Un{X^---XK), (32) 

where K is the number of state changes. 

Lecomte and Tailleur^^-* introduced another method called thermodynamic integra- 
tion. Because the generating function is written as 

/.(A) = iln(e«^*), (33) 

its derivative is 

/^'(A) = ^ = (34) 
which can be interpreted as the average value of q with respect to the population of 
clones. Thus, obtaining numerically and integrating it from to A, we obtain 

/■A 

//(A) = / dA'^^'. (35) 
Jo 

This thermodynamic integration gives a smoother result than the direct evaluation (32). 
In this study, we use both methods and adopt a better result. 

2.5 Known results and the symmetry relations 

There are some known results about the current large deviation for the SEP with 
open boundaries. 
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The cumulants of the current in the SSEP were calculated^^^ under the conjecture 
of the additivity principle as 

(Qb) = J (36) 
{d)c = H (37) 
(..^)c = ^^M^ (38) 

Wb/c - j5 i-^yj 

dpD{p)a{pr-\ (40) 

pb 

Here, the macroscopic parameters are D{p) — 1 and a{p) — 2p(l — p), and the boundary 
parameters are pa — and pb — . Note that all the cumulants are proportional 
to 1/L. 

Next, a fluctuation-theorem-like symmetry was found for the current large deviation 
function for the SSEP3.24) ^ 

fiqn) - fi-qn) = AbQe, (41) 
where As is the constant written as 

.l, = ln^. (42) 

De Gier and Essler^^^ generalized this to the ASEP, where the same equation (41) holds 
with the coefficient Ab given as 

^B = ln^ + (L-l)ln^. (43) 

The coefficient As was derived from the detailed balance condition given by Enaud and 
Derrida,^^^ which is of the form 

In the same manner, we can further extend the symmetry relation for the mean total 
current qa as 

fiqA) - f{-qA) = AaQa, (45) 

where 

Aa = j;^Ab. (46) 
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Fig. 1. Generating function calculated by the largest eigenvalue method when L = 10 and 

those calculated by the Monte Carlo method when L — 50 and 100. The line is interpolated to bring 
clearness to the graph. 

We briefly review other studies related to our study. The WASEP has been well stud- 
ied. Bodineau and Derrida^"^-* assumed the macroscopic fluctuation theory and extended 
the results on the WASEP to the strong asymmetry limit, where they estimated the 
boundary effect on the local Gaussian fluctuation. Depken and Stinchcombe'^^-' calcu- 
lated the joint probability of the density and the current in the TASEP when 7 = 5 = 0. 
Prolhac and Mallick^^-* calculated the cumulants of the current in the WASEP with pe- 
riodic boundary conditions using the Bethe ansatz. 

Now, we see that the large deviation for the total current in the SEP with open 
boundary conditions has not been well studied. 

3. Simulation Results 

3.1 Open boundary SSEP 

First, we show the results on the open boundary SSEP where the hopping rates 
are set a.s pi = pr = l-O- Figure 1 shows some examples of the generating function 
of qA calculated by the two methods; the largest eigenvalue method (L = 10) and 
the Monte Carlo method {L = 50 and 100). The boundary parameters are chosen 
as (a,/3,7,5) = (0.1,0.2,0.9,0.8) just for illustration. The generating functions depict 
parabola-like curves, whose width decreases as the system size increases. 

We fit /iA(A) to the power series up to the sixth order using the least-squares method 
and examine the system size dependence of the coefficients, which are the cumulants of 
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Fig. 2. Comparison between results obtained by the Monte Carlo simulation with L = 50 and fitting 
function of the 6th-order. The broken line shows the quadratic part of the 6th-order fitting function. 

the current. The fitting is performed in the range — 2 < A < 2. The first-order cumulant 
is fixed to the steady current. As seen from Fig. 2, the fitting is successfully performed 
in the entire range with L — 50, whereas the fitting up to the second order significantly 
deviates from the simulation results in the range |A| > 0.5. 

The system size dependence of some low-order cumulants is shown in Fig. 3, where 
the data of 7 < L < 11 arc obtained from the largest eigenvalues and those of the larger 
sizes from the Monte Carlo simulations. The boundary parameters arc (q;,/3,7,5) — 



(0.1,0.1,0.1,0.1), (0.1,0.1,0.9,0.9), (0.1,0.3,0.8,0.5), and (0.3,0.3,0.7,0.7). In every 



case, the cumulants of even order are proportional to the system size but those of odd 
order are not. Thus, the large deviation function cannot be estimated by the mean and 
the variance only; its characterization needs higher-order statistical quantities. However, 
we must note that the present fitting is not very reliable, because the result depends on 
the range of data and on the function to fit with. To obtain more accurate results for 
the higher-order cumulants, we need a wider range of data. In particular, the cumulants 
of odd order have relatively small magnitudes and fitting errors arc significantly large. 
Thus, we could not determine how these cumulants really change with the system size. 
From our results, we conjecture that the generating function split into three parts 

as 



The first term in the rhs represents the mean current, which does not depend on the 




(47) 
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Fig. 3. 2nd- through 6th-order cumulants predicted from the fitting of the generating functions 
given by the simulations. The data of 7 < L < 11 are given by the largest eigenvalue and the other 
data are given by the Monte Carlo methods. The data are obtained using the parameters (a, /3, 7, S) = 
(0.1,0.1,0.1,0.1), (0.1,0.1,0.9,0.9), (0.1,0.3,0.8,0.5), and (0.3,0.3,0.7,0.7). 



system size for the SSEP and is proportional to the system size for the ASEP. The 
second term is the even part, which is proportional to the system size L. We are not 
certain how the odd part behaves. 

We show the data of + ^))/(-^ + 1) = ■j^fJ'A,cvcn{^) when L = 8, 10, 20, 

30, 40, 50, and 70 in Fig. 4, which shows data collapse. This supports the conjecture 
that the even part of the generating function is proportional to the system size. 
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Fig. 4. The function (/i^(A) + /UA(— A))/(Iy + l) is plotted with the system sizes of L = 8(x), 10(+), 
20(*), 30(n), 40(0), 50(A), and 60(v). 



3.2 ASEP with the open boundary conditions 

Now, we present the results on the open boundary ASEP, where the transition 
rates satisfy pi < Pr — 1, and observe whether there are any differences among phases. 
Because the high-density phase is equivalent to the low-density phase due to the particle- 
hole symmetry, we examine the remaining two phases (the low-density phase and the 
maximum-current phase) and the coexisting state of the low- and high-density phases. 
Figure 5 illustrates the generating function /iA(A) in each phase of the ASEP with 
system size L = 10 and pi = 0.5. We note that no corresponding analytical or numerical 
results are indicated in the literature. It is remarkable that any qualitative differences 
are not observed among the three states. As in the SSEP, we have attempted to fit 
the data with the 6th-order polynomial, and the result for L — 50, pi = 0.5, and 
(q;,/3,7, (5) = (0.9,0.9,0.1,0.1) is shown in Fig. 6. The fitting seems to work well if a 
certain neighborhood of the minimum is excluded. We confirmed this property also in 
other phases. This fiattening of the minimum has been reported also in a system with 
periodic boundary condition. Clearly, the generating function cannot be assumed to 
be quadratic. 

Similarly to the case of the SSEP, we may divide the generating function for the 
ASEP into three parts as in eq. (47). We illustrate (/iA(A) + /iA(— A))/(L + 1) in Fig. 
7 for various system sizes L = 8, 10, 20, 30, and 40. In relatively small systems, the 
maximum current phase shows very good data collapse, while the degree of data collapse 
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Fig. 5. Generating function /xa(A) in the ASEP when L = 10 and pi = 0.5. Crosses, pluses, and 
stars correspond to the coexisting phase ((a, /3, 7, 6) = (0.1, 0.1, 0.1, 0.1)), the maximum-current phase 
((0.9,0.9,0.1,0.1)), and the low-density phase ((0.1,0.9,0.9,0.1)), respectively. 




Fig. 6. Fitting of the generating function /za(A) with 6th-order function in the maximum-current 
phase. We may fit /Ua(A) well in the region outside the minimum, and this can be seen in other phases. 



is not SO good in other phases. In every case, however, the data seem to converge in 
the large-system hmit. Thus, as in the SSEP, the even part of the generating function 
is proportional to the system size when it is sufficiently large. 

The large deviation function is obtained from the generating function via the Legen- 
dre transformation. Figure 8 shows the plot for L = 10, {a,/3,j,6) = (0.1,0.9,0.9,0.1), 
and pi = 0.1, 0.2, 0.3, 0.5, all of which correspond to the low density phase, in the range 
—4 < q < 4:. A blowup near the minimum of f{qA) is shown on left of Fig. 8. We 
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Fig. 7. The function {haW + A*a(— A))/(L + 1) is plotted with the system sizes of L = 8(x), 
10(+), 20(*), 30(n), and 40(0)) in the low-density phase (top), maximum-current phase (middle) and 
coexisting phase (bottom). 
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Fig. 8. Large deviation funetion given by tlie Legendre transform in tlie range —5 < cja l£ ^- Tlie 
parameters are L = 10, {a,(3,j,S) = (0.1,0.9,0.9,0.1) (low-density phase), and p; = 0.1, 0.2, 0.3, and 
0.5. The derivative of the left figure is shown in the right figure. 



note that a cusp appears near g = and its sharpness increases as pi decreases. This is 
confirmed by the plot of the derivative of the large deviation function shown on right of 
Fig. 8. The cusp means that the probability of positive current rapidly decreases with 
L. Such a cusp is seen also in other phases. In the limit — > 0, the particles do not 
move leftward except at the boundaries. Thus, the large deviation function is expected 
to diverge for > in the limit. This is why the cusp is generated. 



3.3 Comparing the two currents 

Thus far, we have presented the simulation results for the total currents. In this 
subsection, we compare the characteristic functions of and qs- First, we introduce 
the dimensionless current qx defined by 

where (qx) is the expectation, and then we write the generating function iJ>xi^) of the 
current qx- We compare A*^(A) and A*s(A) in the simulations. 

The dimensionless current defined here is effective when (qx) is known. Fortunately, 
the steady current is exactly calculated for both the SSEP*^®) and the ASEP^") as 

5 a 

- ^ (AQ\ 
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Fig. 9. Plot of /U^(A)(x) and Ms (A)(+) for three sets of parameters. We observe a clear coincidence 
of /x^(A) and Ms (A). 



for the SSEP, and 

^ (p; — l)po(l — Pa), low-density and coexisting phases 
{<Ib)asep = ipi — — Pb), high-density phase 

y maximum-current phase 

for the ASEP with pi < = I, where pa and pb are given as 
1 



Pa 



and 



1 + a 
b 



Pb = 



a — 



6 = 



l-pi-a + 'y + - pi - a + 'yy + Aa^ 



2a 



l-pi-5 + (5 + ^{l-pi-5 + (5f + AI35 



(50) 



(51) 



(52) 



1 + 6 2^ 
Henceforth, we omit the suffix 'SSEP' or 'ASEP' from the steady current (gx), for it 
can be easily understood from the context. 

The simulation results from the method of largest eigenvalues are shown in Fig. 9. 
It is clear that /U^(A) and p^{\) coincide with each other. 

However, the two large deviation functions show differences in the Monte Carlo 
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Fig. 10. Plot of n^{X) and /i^(A) for the SSEP and the ASEP. The inset shows the extended plot 
of the coinciding area. 



simulation, as shown in Fig. 10. They agree with each other only when |A| is very 
small and differ greatly outside the region. This difference means that /i^(A) is not 
appropriately calculated in the Monte Carlo simulation. One reason for this is the lack 
of statistics. Because qb counts only the flux at the left boundary, the number of such 
events is much smaller than that in the case of qa- Another possibility is that the 
number of clones is not sufficiently large. However, no signiflcant change has been seen 
in a simulation where we doubled the number of clones and the time duration, as seen 
in Fig. 11. 

Thus, the difference between the two currents can affect results of the Monte Carlo 
simulations, though they represent essentially the same physical quantity. The use of 
the total current qa has an advantage in this respect. 
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Fig. 11. Plot of ij-bW for different time step K and clone number A'' values. 



3.4 Symmetry relations 
The symmetry relation 

f{qA) - fi-QA) = Aqa (53) 

should hold in the SEP. First, we show the results for small system sizes obtained from 
the direct evaluation of the largest eigenvalue. We plot /{qa) — /{—Qa) and AaQa with 
Aa given by eqs. (43) and (46) in Fig. 12. The parameter sets used are (a, ^,7, 5) = 
(0.1,0.9,0.8,0.2) for the low-density phase that we denote by x, (0.1,0.1,0.1,0.1) for 
the coexisting phase that we denote by +, and (0.9,0.9,0.1,0.1) for the maximum- 
current phase that we denote by *, and the other parameters L = 10, pi = 0.5, and 
Pr = 1.0 used are the same. The solid line shows AaQa- We see that the symmetric 
relation is well satisfied in this case. 

Next, we consider the results of the Monte Carlo method. In this case, we used 
the thermodynamic integral to compute the generating function and obtained large 
deviation functions via the Legendre transform. In Fig. 13 we plot /(^a) — f{—QA) for 
the SSEP and the ASEP. The parameter sets are given in the inset of the figure. A 
neghgible deviation is observed for the SSEP, while a small deviation is apparent in the 
ASEP. 

The deviation increases with the system size. The error may stem from the insuffi- 
cient number of clones or the resolution of the data, which affects the precision of the 
Legendre transform. For large system sizes, the number of necessary clones becomes so 
large that we could not achieve a higher accuracy. Thus, deviations in Fig. 13 should 
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Fig. 12. Plot of /{qa) — /{—Qa) from the direct evaluation results represented by x, +, and *, and 
that of AaQa from the analytical result obtained using eq. (46), which is shown by the solid line. 

be interpreted as representing the precision and limitation of the present Monte Carlo 
simulations. 

3.5 Convergence problem in the Monte Carlo method 

The Monte Carlo method sometimes presents non-negligible errors. It happens when 
the hopping asymmetry and the system size are large. For example, Fig. 14 illus- 
trates the instability observed in the numerically obtained generating function. Here, 
we have chosen the parameters L = 50, (a,/3,7,5) = (0.9,0.9,0.1,0.1) (maximum- 
current phase), and pi = 0.2. This figure shows characteristic fiuctuations in the region 
0.2 < A < 1.5. The true generating function should be a smooth or at least convex 
function. However, the obtained data docs not satisfy the expected convexity. To find 
the cause of the fiuctuations, we investigated the initial condition dependence of our 
Monte Carlo simulation. In Fig. 15, we show how t~^ ln(Xi • • • Xk) evolves with time in 
the cases of A = — 1 and A = 1 each with three different initial configurations, namely, 
the empty initial condition (all = 0), the random initial condition (tj is randomly 
chosen with probability 1/2), and the half-filled initial condition (tj = 1 for i < L/2 
and Ti — for the other value oi i). The other parameters used are the same as those 
indicated in Fig. 14. Note that A = — 1 corresponds to the stable region and A = 1 the 
unstable region. The top of Fig. 15 shows the case A = — 1, where the three samples 
exhibit convergence to a common value at t = 10"^. The bottom shows the case A = 1, 
where the convergence is very slow and the initial-condition dependence remains even 



J. Phys. Soc. Jpn. 



FULL PAPERS 





a= 


=0.1,P=0-1.'F0.9,6=0.9,L=50 + 


0.4 




=0.1,P=0-1.'F0-9,6=0.9,L=30 x . 






(-log(3)*4)/51*x — 






/ lrt/i/0\*/l\/Oi *v 

(-iog(oj 4j/oi X — 





f(q.)-fK) 


SSEP 








0.4 













-4 -2 2 4 



4 




-4 -2 2 4 



Fig. 13. Comparison of /{qa) — .f{—lA) values for tlic SSEP (top) and ASEP (tire others). For the 
SSEP, the system sizes of L = 30 and 50 are plotted. For the ASEP, the plot of L = 30 is given in the 
second figure and the plot of L = 50 is given in the third figure. The parameters are given in the inset. 
The solid and dotted lines give the analytical results. 
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Fig. 14. /x(A) when L = 50, {a, (3,^,5) = (0.9,0.9,0.1,0.1) (maximum-current phase), and pi = 0.2. 

at t = 10^. This slow convergence is considered to be related to the instabihty shown 
in Fig. 14. 

The convergence property should be reflected in the spectrum of eigenvalues of the 
A- modified transition matrix W'^. Thus, we calculated the second largest eigenvalue ^2- 
In Fig. 16, we show the system size dependence of the gap ^max — C2 for several A values 
in the systems of sizes L — 6, - ■ • ,13, where Cmax is the largest eigenvalue. The top figure 
shows the cases A = —0.1, 0.0, 0.4, and 1.6 on the log-log scale and the bottom figure 
shows the cases A = 0.0, 0.1, 0.16, 0.2, and 0.24 on the semi-log scale. The figure shows 
that the magnitude of the gap is proportional to a certain negative power of L in the 
stable region, whereas the gap decreases faster than the exponential decrease with L in 
the unstable region. In the bottom of Fig. 16, we sec that the system size dependence 
of Cmax — C2 changes from the power to exponential and finally to stretched exponential- 
like. Thus, one reason for the instability of the Monte Carlo calculation found here is 
considered to be that ^2 becomes close to Cmax as L increases. This behavior suggests 
that a certain type of phase transition may occur as reported recently.^^'^^^ 

4. Discussion and conclusions 

Our simulation results show A*^(A) — )U^(A), which implies 



From the calculation given by Bodineau and Derrida for the SSEP,^^) we may assume 
{Qb)c ~ X' tii-ViS (?2l)c ~ L"'~^. This relation does not coincide with our result 




(54) 



J. Phys. Soc. Jpn. 



FULL PAPERS 



20 
^(-1) 

10 



1 10 100 1000 10000 100000 

Time step 

-1 

mi) 

-2 
-3 
-4 



1 10 100 1000 10000 100000 

Time step 

Fig. 15. Time evolutions of m(— 1) (top) and (bottom) when L = 50 andp; = 0.2 starting from 
three initial conditions: empty, random, and half-filled. 



{q'a')c ~ L. These conflicting results suggest the necessity of further study. The calcu- 
lation by Bodineau and Derrida assumes the additivity principle, the scaling relation 
IXqe) ~ 9i(lBL)/L, and the large deviation locally defined to be quadratic. Our results 
do not support the scaling relation and the Gaussian property of the local large devi- 
ation function. Both assumptions may need revision. However, our results have some 
drawbacks. The estimation by Bodineau and Derrida is given in the entire parameter 
region, while our simulation is limited to several sets of parameters. Furthermore, the 
present Monte Carlo method has errors in several cases. The improvement of the Monte 
Carlo method to obtain the generating function is necessary. We must note that the 
DMRG method is another method for obtaining a large deviation,^^) and not only the 
continuous time, which is well studied, but also the discrete time case based on the 
DMRG method for the steady state^^^ can be established. 

The relation qa ~ Lqb can be derived from the plausible argument considering 
the time-dependent and time-independent contributions to the current. The details are 
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Fig. 16. Cmax — C2 versus the system size L = 6, - ■ ■ ,13. Tlie top figure sliows tlie cases A = —0.1, 
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presented in § 2. Actually, the two generating functions for the two currents coincide 
with each other in relatively small systems, as shown in Fig. 9. Although they are 
considered to be the same also in large systems, we could not confirm it because the 
Monte Carlo calculation is not suitable for determining qb- Thus, in a practical sense, 
they have differences. 

We observed the cusp around q — in the current large deviation function in 
the ASEP. A similar cusp in the large deviation function for entropy production was 
reported by Mehl et al.^^^ in a system of an overdamped Brownian particle under the 
periodic potential and uniform external field. They attributed the generation of the 
cusp to the sublinear response of the particle to the external field, which was derived 
from the potential within the site. In our study for the ASEP, the subhnear response 
may be caused by the interactions of particles, where the cause is different from that 
in the case of a one-particle system. The physical basis of the cusp we consider is that 
the probability of observing the current toward a lower hopping rate becomes 0, which 
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means that the current becomes zero in that direction. This explains why the cusp is 
observed, though this does not explain how the cusp is generated. 

In conclusion, in this study we have calculated numerically the current large devi- 
ation function in the SSEP and ASEP with open boundaries. We have found that the 
generating function and the large deviation function deviate from the quadratic form. 
The system size dependence of the even order cumulants is found to be proportional 
to L in both the SSEP and the ASEP. The generating functions defined by different 
definitions of the current appear to be the same for small systems. For the large system 
size, the Monte Carlo simulation gives unsatisfying results for the current through a 
boundary, from which we conjecture that the use of the total current is suitable for the 
Monte Carlo calculation of the large deviation function. The symmetry relation for the 
total current is calculated both numerically and analytically, and the obtained results 
are in good agreement. For the Monte Carlo calculations, a small deviation is observed 
in the symmetry relations and the deviation depends on the boundary parameters. We 
found that the convergence of the Monte Carlo simulation for calculating the large de- 
viation function in the ASEP depends on the system size and the hopping parameter 
of a particle. The convergence is not good when the system size is large and the asym- 
metry of the hopping is large. In the region where the convergence is not good, the 
difference between the largest and second largest eigenvalues of the modified transition 
matrix depends on the system size as the stretched exponential, which is faster than 
the exponential decay. 
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